Spatial, temporal and molecular dynamics of swine influenza virus-specific CD8 tissue resident memory T cells

For the first time we have defined naïve, central memory, effector memory and differentiated effector porcine CD8 T cells and analyzed their distribution in lymphoid and respiratory tissues after influenza infection or immunization, using peptide-MHC tetramers of three influenza nucleoprotein (NP) epitopes. The hierarchy of response to the three epitopes changes during the response in different tissues. Most NP-specific CD8 T cells in broncho-alveolar lavage (BAL) and lung are tissue resident memory cells (TRM) that express CD69 and downregulate CD45RA and CCR7. NP-specific cells isolated from BAL express genes characteristic of TRM, but gene expression differs at 7, 21 and 63 days post infection. In all tissues the frequency of NP-specific CD8 cells declines over 63 days almost to background levels but is best maintained in BAL. The kinetic of influenza specific memory CD8 T cell in this natural host species differs from that in small animal models.


INTRODUCTION
Immunity to influenza A viruses (IAV) has been intensively studied over many years and the role of neutralizing antibody in protection against homologous virus is well established (reviewed 1 ). Moreover experimental studies in mice and humans have revealed the importance of pre-existing cellular immunity in heterosubtypic protection (reviewed 2 ). The development of live attenuated influenza virus vaccines administered by nasal spray, now commercially available for humans and pigs, followed from the discovery that many lymphocytes reside in non-lymphoid tissues and that local immune responses are important in protective immunity 3,4 .
Intravenous infusion of anti-lymphocyte antibodies prior to isolation of lymphocytes from tissues has led to the definition of populations of lymphocytes that are considered to be tissue resident because they do not enter the bloodstream and therefore are not stained by the infused antibody. The majority of these cells have a memory phenotype and are therefore termed tissue resident memory cells (TRM) [5][6][7][8][9][10] . In mice TRM have been shown to exceed the number of T cells in the lymphoid system and to play important roles in maintaining local immune memory 11 . Adoptively transferred lung TRM cells from mice that were previously infected with Sendai virus conferred protection to naïve animals, demonstrating the essential contribution of this population to immunity 12 . TRM are the predominant population in the adult human lung and antigenspecific cells were found at stable frequencies years after pathogen encounter, indicating their key role in respiratory infections 13 . Mouse and human TRM express CD69 and more variably CD103, both facilitate tissue retention 6,[14][15][16][17] . CD103 expression on lung TRM is dependent on TGFβ signalling 18 .
Porcine physiology closely resembles that of humans, pigs have a similar distribution of sialic acid in their respiratory tract and are infected with similar influenza viruses, making them a powerful large animal natural host model to study immunity to influenza [19][20][21] . Reassortment of influenza gene segments in pigs is well established, making pigs an animal reservoir of influenza viruses that pose a threat of zoonotic infections in humans 22 . Influenza is also a significant economic burden to farmers 22 . We have shown that pigs and humans make very similar antibody responses to influenza viruses and that local respiratory tract immune responses play an important role in protective immunity. Following immunization with influenza vaccines or influenza infection, there are high frequencies of influenza-specific T cells in broncho-alveolar lavage (BAL), lung tissue and tracheobronchial lymph nodes (TBLN) [23][24][25][26] . However, the lack of monoclonal antibodies (mAbs) to porcine CD69 and CD103 has made studies of porcine TRM challenging, although we have previously identified porcine TRM, following infusion of anti-CD3 antibody 25,27 . Nevertheless, their phenotype and function during influenza infection remain poorly defined and there have been few ex vivo studies of antigen-specific T cells. In addition, while the phenotype of porcine helper T cells has been thoroughly analyzed [28][29][30][31] , CD8 T cells are less well characterized and CCR7, essential for migration to lymphoid organs and CD45RA, a known differentiation marker, both widely used in human immunology, have never been studied in combination to define subsets of porcine T cells.
Here we examined IAV-specific CD8 TRM throughout the respiratory tract. We focussed on CD8 TRM in BAL as this population of airway cells is relatively accessible in many species including humans, contains almost exclusively TRM apart from alveolar macrophages, and does not require extraction procedures that might alter the cellular composition and phenotype. Here, for the first time we showed that CD45RA and CCR7 together identify porcine CD8 T cell subsets similar to those in humans and we described the expression of CD69 in CD8 T cells in several tissues, using a newly generated antibody 32 . We examined antigen-specific CD8 TRM in the context of influenza infection and immunization in inbred Babraham pigs 33 , using peptide-SLA tetramers carrying a previously identified 34 and two novel influenza nucleoprotein (NP) epitopes. This enabled us to define the phenotype of influenzaspecific CD8 T cells, analyze their distribution in different tissues, define their transcriptional profile at different times after infection and model the dynamics of the response.

Subsets of porcine T cells and identification of CD8 TRM
Antibodies to CD45 isoforms and CCR7, which distinguish T cell subsets in several species, were used for initial characterization of CD8 T cells from different tissues of Babraham pigs. BAL, which contains~80% macrophages as in other species, is of particular interest as it contains airway T cells that are at the frontline of protection against respiratory pathogens 35 . Among the lymphocytes in BAL~20% were CD4 T cells,~20% CD8 and a slightly higher proportion γδ T cells 24 .
Because porcine T helper cells express CD8α when activated 30,36 , we used an antibody to CD8β combined with CD45RA and CCR7 antibodies to identify and characterize CD8 T cells. Four populations were apparent (Fig. 1a) in blood, spleen and TBLN, which appear to correspond to those defined in humans as naïve (CD45RA + CCR7 + ), central memory T cells (TCM) (CD45RA -CCR7 + ), effector memory T cells (TEM) (CD45RA -CCR7 -) and differentiated effector T cells (TDE) (CD45RA + CCR7) 37 . The mucosal tissues (lung and BAL) contained higher proportions of cells with TEM and TDE phenotypes. To confirm this differentiation scheme for CD8 T cells, we examined expression of CD27, a marker of less differentiated cells and perforin, expressed by effector CD8 T cells. These experiments indicated that the proposed naïve and TCM cells in all tissues expressed CD27 but little or no perforin, while TEM and TDE cells were heterogeneous, expressing one or the other marker or neither, as in humans (Supplementary Fig. 1) [37][38][39] . These data support the identification of naïve and TCM cells and indicated that TEM and TDE are more differentiated cells with effector function. To confirm this, we examined PMA/ionomycin induced cytokine production of the four subsets of CD8 T cells, sorted from peripheral blood mononuclear cells (PBMC). A high proportion of TEM produced IFNγ and TNF while lower proportions of TDE and TCM and few naïve cells did so (Fig. 1b, c, d). Naïve and TCM cells secreted mainly TNF (4.4% and 8.7% respectively) while many TEM were double producers (IFNγ + TNF + 18.9%). TDE produced predominantly IFNγ (8%) with a smaller proportion secreting both IFNγ and TNF (4.9%).
We have previously shown that BAL and a large proportion of lung T cells are not stained by intravenous anti-CD3 antibody 25,27 , indicating that these populations are TRM. BAL T cells are almost exclusively highly differentiated. Staining with CD69, a marker of activation and tissue residency 40,41 showed that, as in humans, CD69 was absent or minimally expressed on blood CD8 T cells, while the highest expression was found in BAL and TBLN. In TBLN CD69 expression may indicate either recirculation of TRM 10 or recent activation of T cells. CD8 TDE in the lungs expressed lower levels of CD69 (Fig. 1e).
These data demonstrated that in pigs, CD69 expression is tissue dependent on activated/memory CD8 T cells. BAL CD8 T cells are predominantly TEM phenotype with high expression of CD69 and lack of staining by intravenous CD3 antibody, indicating that they are TRM. TEM are heterogenous in expression of CD27 and perforin and produce high levels of cytokines on stimulation.
Kinetics and phenotype of influenza-specific CD8 T cells Having established the phenotype of CD8 T cells from unimmunized pigs, we wished to define the temporal dynamics of the immune response to IAV. Inbred Babraham pigs were infected with H1N1pdm09 in four experiments and pigs culled from 1 to 63 days later 24 (Fig. 2a). We determined the proportions of influenza NP-specific CD8 T cells binding to tetramers carrying the previously defined DFE epitope and tetramers of the newly identified AAV and VAY epitopes ( Supplementary Fig. 2) in various tissues over time (Supplementary Table 1), starting at 6 DPI when a cellular response is first detectable 24 . The results and curves fitted to the data are shown in Fig. 2b and indicated that the peak proportion of CD8 + T cells specific for each tetramer and the timing of the peak differed amongst tissues and between tetramers. The highest response for all tetramers, expressed as a percentage or absolute count, was found in the BAL, followed by lung and local lymph nodes, ( Fig. 2b and Supplementary Fig. 3). Interestingly, the modelled response in PBMC peaked earlier compared to BAL, lung or nasal turbinates (NT) ( Table 1), in accordance with the idea that cells generated in lymph nodes traffic to local tissues via the blood. AAV responses peaked 3 to 17 days later than the peak of DFE or VAY responses in different tissues, while DFE and VAY shared similar kinetics in most tissues except PBMC, where DFE peaked earlier (Table 1 and Fig. 2b). The magnitude of AAV responses was greater than those to VAY and DFE in all tissues, while DFE was higher than VAY in all tissues except for TBLN (Table 1 and Fig. 2b). At the later time points of 21, 42 and 63 DPI we also examined tracheal CD8 T cells and here too the AAV response was higher than that to DFE and VAY (p = 0.03 at 42 DPI), with DFE being significantly higher than VAY (p = 0.03) (Supplementary Fig. 4A, B).
We next analyzed changes in the proportions of the different tetramer + cells in each tissue by fitting curves to the data (Supplementary Table 2). BAL, lung and spleen showed similar changes in proportions of the different tetramer + populations. Initially AAV was lower than VAY and DFE in all tissues but by 30 DPI AAV was dominant (60%), while in most tissues DFE and VAY declined (~20%) (Fig. 3a). However, in NT the frequency of cells specific for DFE remained constant (25%), while that for VAY decreased and that for AAV increased. In PBMC the proportions changed only gradually (Fig. 3a), though the small number of PBMC data points indicate that this observation should be interpreted with caution.
Influenza specific T cells in the nasal mucosa of mice have been shown to be longer lived and decline less rapidly than those present in the harsh environment of the lung 42 . Using the curves fitted to the proportion of different tetramer + cells in each tissue (Fig. 2b), we investigated the decline of tetramer + T cells in all tissues after the peak response in more detail (Fig. 3b, c). The decay in the proportion of cells specific for AAV was slower than that for cells specific for either DFE or VAY in all tissues. Furthermore, the proportion of cells specific for AAV decayed most slowly in BAL.
Collectively, our results shows that the frequency of different tetramer + T cells varies between tissues, with the highest frequency in BAL. AAV tetramer + cells are dominant at later time points. In general, responses in PBMC peaked earlier compared to local tissues but waned more rapidly and did not reflect events in mucosal tissues. We did not observe a more rapid decay in lung and BAL compared to NT, as previously reported in mice 42 .
BAL T cells maintain a stable surface phenotype but transcription alters over time We next analyzed the phenotype of the NP-specific T cells present locally (in BAL and TBLN) and systemically (PBMC) (Fig. 4a). The majority of tetramer + cells were TCM or TEM throughout the time course in all tissues analyzed. There was a slow but steady increase in the proportion of TCM with time in PBMC (Fig. 4b). Most influenza-specific T cells in local lymph nodes were TCM (80%) (Fig. 4a, b). In contrast, in BAL the majority of cells are TEM (78.5%), with only a small number of TCM (20.3% at 21 DPI) (Fig. 4a, b). Changes in TDE and naïve cells are not shown as the numbers of these cells were too small for reliable analysis.
We used RNA sequencing (RNA-seq) to compare the transcriptome of CD8 T cells specific for the previously defined DFE NP epitope 34 . DFE + T cells were isolated from BAL by cell sorting at different time points (7, 21 and 63 DPI) (Fig. 5a). Differential gene expression analysis was applied to compare these three groups. 4,666 genes were expressed at significantly different levels (p adj value ≤ 0.05 and |log 2 fold change | > 1) at 7 DPI versus (vs) 21 DPI while 1,198 were differentially expressed in 7 DPI vs 63 DPI and only 560 in 21 DPI vs 63 DPI (Fig. 5b). At 7 DPI several upregulated genes were involved in cell growth, movement (Igfbp2) and proliferation (Ctla4,Kif11, Kif18a, Shmt1) while at 21 DPI genes linked with T cell activation (Tagap, IL2ra, Csf2, Dgkg) and adhesion (L1cam, Cass4) were highly expressed (Fig. 5c). Interestingly, comparison of 7 DPI and 63 DPI revealed differential expression of genes involved in metabolism (Jazf1, Atp8b4, Igf2bp3,), transcription factors (Litaf), T cell development and proliferation (Ccnd3,   Table 3). The TGFβ signalling pathway, known to be involved in mucosal residency, was also upregulated at 21 DPI.
We examined key gene expression features of TRM, previously identified in the human lung 13,43,44 (Supplementary Fig. 5). BAL cells from 63 DPI upregulated a gene related to integrins (Itga1), the TRM transcription regulator gene Znf683 and downregulated genes involved in migration (Sell, S1pr1) as in humans. In addition, at the earlier timepoint of 7 DPI DFE + T cells expressed more cytotoxicity related genes (GzmA, GzmH, Prf1 and Ccl5) while starting from 21 DPI genes involved in cytokine signalling and secretion were upregulated (Ifng, Tnf, Il13, Tgfb1 and Tnfsf13b) ( Supplementary Fig. 5). Interestingly, CD69 gene expression changed with time, peaking at 21 DPI, as did CD103 gene expression (Itgae) with a peak at 7 DPI.
Despite the similar phenotype, these data suggests that gene expression in CD8 T cells at the site of infection changes over time, with genes involved in proliferation and migration upregulated at 7 DPI while cytokine related pathways are upregulated at 21 DPI.
BAL tetramer + cells are a highly differentiated population As RNA-seq analysis revealed changes in CD69 gene expression with time, we analyzed CD69 protein expression in BAL, TBLN and PBMC by fitting curves to the data (Supplementary Table 4 and Fig. 6a). CD69 expression decayed only slightly with time in BAL for all tetramers and also for AAV and VAY labelled T cells in TBLN, while it decayed faster in TBLN DFE + T cells (0.022/day for DFE, 0.003/day VAY and 0.001/day AAV) ( Table 2 and Fig. 6a).
Pathway analyses revealed upregulation of genes related to Th1/ Th2 differentiation of CD4 T cells in BAL CD8 T cells at 21 DPI compared to 7 DPI. We validated these findings by analyzing and modelling the expression of the transcription factor T-bet, involved in Th1 differentiation and homing to inflammatory sites, and Eomesodermin (Eomes), involved in induction of memory and effector T-cell differentiation 45 (Fig. 6b). The highest expression of this transcription factor was found in AAV + T cells in TBLN followed by VAY and DFE with little change over time, while similar expression was present in PBMC for all tetramers, with no decay ( Fig. 6b and Table 2). Eomes was poorly expressed in BAL, only detectable at early time points and decayed rapidly in all tetramer + cells ( Table 2). In contrast, T-bet expression differed greatly among tissues and tetramers. T-bet was highly expressed in TBLN and PBMC, and decayed more slowly in TBLN for all tetramers. In TBLN, T-bet decreased more gradually in AAV than DFE or VAY responding T cells ( Fig. 6b and Table 2). There was low expression of T-bet at early times in BAL, which declined rapidly and was undetectable at later time points in all tetramers (Table 2).
These data suggested that BAL TRM may have already switched off Eomes and T-bet protein expression and are no longer undergoing active Th1/Th2 differentiation. We therefore analysed Ki67 expression as a proxy for cell proliferation which is normally linked to differentiation. Whereas high frequencies of Ki67 + tetramer binding cells were found in PBMC and TBLN at early time points, only 14% of BAL cells were Ki67 + at 6 DPI and Ki67 expression was barely detectable at 63 DPI (1.8%) (Fig. 6c).
The RNA-seq data indicates changes in expression of many genes related to cytokine production over time. Lymphocytes isolated from PBMC, TBLN and BAL were therefore stimulated with H1N1pdm09 and the production of IFNγ, TNF and IL-2 by tetramer binding cells was assayed by intracellular cytokine staining. We compared the responses of the dominant responding AAV population with DFE cells which decline more rapidly. Despite high expression of T-bet, PBMC tetramer + cells produced a limited amount of IFNγ (5.9% in DFE + , 4,3% AAV + at 7 DPI) which was almost undetectable after 21 DPI (Fig. 7a, b). Similar kinetics were found in TBLN, with IFNγ and IFNγ/TNF co-producing cells being the most abundant. The highest responses were in BAL, with consistent production of cytokines (predominately IFNγ and TNF) even at 63 DPI. No significant difference was observed between DFE + and AAV + populations, but there was a trend toward a higher proportion of triple producers (INFγ + TNF + IL-2 + ) at later time points in AAV + cells compared to DFE + (Fig. 7b).
In conclusion, despite changes over time in transcription of CD69 and T cell differentiation genes, we did not find corresponding differences in the level of protein expression of CD69, T-bet and Eomes in BAL. In contrast TBLN and PBMC expressed T-bet and Ki67 at early time points and Eomes up to 63 DPI. Ex vivo stimulation with H1N1pdm09 resulted in IFNγ and TNF cytokine secretion at 7 DPI, while upregulation of related genes was highest at 21 DPI. Irrespective of changes in gene expression, BAL TRM maintained a stable surface phenotype and produced abundant cytokines for at least 63 DPI.
Aerosol immunization generates a powerful CD8 response with a similar phenotype to influenza infection To investigate whether immunization elicits similar responses to natural infection we administered S-FLU, a single cycle influenza vaccine, by aerosol to Babraham pigs and boosted them 3 weeks later. In a first experiment three pigs were culled 3 weeks post boost (WPB) while in a second, S-FLU immunized pigs (n = 6) were challenged with H1N1pdm09 and culled four days later (Fig. 8a). Control pigs were left untreated in the first or were challenged without prior immunization in the second experiment. Anti-porcine CD3 mAb was administered intravenously (i.v.) 10 min prior to sacrifice to distinguish between tissue resident (CD3 i.v. -) and circulating T cells (CD3 i.v. + ) ( Fig. 8a and Supplementary Fig. 6). We enumerated S-FLU specific cells using the DFE tetramer, as this epitope is conserved in the vaccine and challenge viruses. In BAL, 22.6% of CD8 T cells were DFE + at 3 WPB but only 5% 4 days post challenge (DPC) (Fig. 8b). Conversely, there were higher numbers of TBLN and PBMC DFE + cells 4 DPC compared to 3 WPB (0.8% vs 0.3% in TBLN and 0.2% vs 0.1% in PBMC) (Fig. 8b). As we have shown before, BAL and TBLN cells were completely inaccessible to blood (CD3 i.v. -) while most spleen cells were labelled with the infused antibody ( Supplementary Fig. 6), confirming that BAL cells are truly tissue resident and that lymph node cells are also outside the blood stream 25,27 .
We then used CD45RA and CCR7 to study differences in phenotype. Three weeks after S-FLU immunization, BAL cells were almost exclusively TEM as after natural infection and this did not change after challenge (Fig. 8c). Following challenge however, there was a rapid increase in the proportion of TCM in TBLN and PBMC (from 33% at 3 WPB to 77% at 4 DPC in TBLN and 0% to 46% in PBMC).
In the earlier H1N1pdm09 infection experiment, relatively few BAL CD8 cells expressed Ki67. We therefore examined whether DFE + T cells expressed Ki67 early after immunization and challenge with live virus. Ki67 was minimally expressed at 3 WPB in all tissues, and at 4 DPC in BAL. In contrast a high proportion of DFE + TBLN cells and DFE + PBMC expressed Ki67 after challenge (Fig. 8d). As during influenza infection, 3 weeks after the S-FLU boost, DFE + BAL cells lacked T-bet and Eomes but Eomes was expressed at a low level in PBMC, and T-bet minimally in TBLN and in a low proportion of PBMC. At 4 DPC, the distribution of Eomes was unchanged but T-bet expression increased in all tissues including BAL (Fig. 8d).
Taken together these results suggest that S-FLU immunization by aerosol elicited a stronger T cell response in BAL compared to influenza infection but the responding T cells had a similar differentiated, largely non-proliferating phenotype. The proportion of DFE tetramer + cells in BAL decreased four days after challenge with infectious virus perhaps due to activation induced cell death following encounter with antigen in the epithelium.

DISCUSSION
In this study, we have described in detail for the first time porcine influenza-specific CD8 T cells in lymphoid and non-lymphoid   2). b Decay from the peak of response within the indicated tissue. c Decay of tetramer+ CD8 T cells in different tissues (as indicated). Note: DFE and VAY have the same dynamics in BAL, VAY has the same dynamics in NT and PBMC and AAV has the same dynamics in NT and spleen so these decay curves overlap.
other species will require investigation of their proliferative potential, telomer length and expression of survival and senescence related genes (reviewed in 46 ). We have already shown that the majority of TBLN, BAL and, to a lesser extent, lung T cells are inaccessible to intravenous anti-CD3 antibody, it was therefore of interest to examine the expression of CD69, a marker of activation and tissue residence, on these cells. CD69 is a lectin that binds the sphingosine 1-phosphate receptor-1 (S1PR-1). S1PR-1 binding limits responsiveness to sphingosine 1phosphate, a mediator of egress to the bloodstream, and therefore promotes tissue retention. As in other species, we found that CD69 is minimally expressed on blood T cells although upregulated if PBMC are activated with a mitogen (data not shown), but in BAL, CD69 is expressed on a high proportion of CD8 T cells. In TBLN and spleen CD69 + T cells were also present, perhaps indicating that a proportion of cells in these organs may have been recently activated. A high proportion of NP specific T cells in TBLN expresses CD69 from 6 days after influenza challenge until 63 DPI (Fig. 6a) perhaps representing either continuing antigenic stimulation or recirculation of TRM 10 . Murine studies identified a population of CD4 and CD8 TRM, inaccessible to the bloodstream, lacking CD69 and CD103 expression, highlighting the heterogeneity of this population 17 . Other markers such as CD49a and the chemokine receptor CXCR6, which promote tissue retention of TRM, should therefore be studied in the future to further characterise porcine TRM 47,48 .
Most studies on T cells responses to IAV in pigs have relied on ex vivo restimulation with live virus, with consequent changes in phenotype 30,49,50 . Here, we reported the the identification of two new CD8 epitopes in influenza NP, which allowed us to compare changes in the distribution, phenotype and gene expression of different antigen specific cells from 6 to 63 days post influenza infection 34 . Intriguingly, the modelled kinetics of the responses to the three epitopes are very different, with DFE tetramer numbers showing a very early peak in PBMC and spleen but peaking in all other tissues at similar times to VAY. The response to the third epitope AAV appears later, but by 20-30 DPI it is dominant in all tissues and persists for longer than the response to the other two epitopes. Others have reported marked differences in the hierarchy of the responses to immunodominant epitopes of influenza in mice 51,52 . Pizzolla and colleagues have shown that T cells specific for immunodominant epitopes in NP, PA and PB proteins of influenza in the nasal mucosa shared a similar hierarchy with systemic responses while in the lung all immunodominant epitopes were equally represented 42 . In our model, despite an initial peak of DFE + CD8 T cells, at the memory stage, the AAV epitope dominates in all tissues, with no marked difference between the organs. Furthermore, the response of tetramer positive cells does not decline more rapidly in the lung or BAL than blood or spleen, although in mice both of these have been postulated to be hostile environments inducing transcriptional and epigenetic changes that promote apoptosis 53 . Further studies are required to understand why AAV responses are immunodominant compared to DFE and VAY. Immunodominance may be linked to peptide-MHC affinity and it has also been reported that early IFNγ production can provide an advantage to a given antigen specific population 54,55 . However, we did not observe significant differences in cytokine production of AAV + and DFE + CD8 T cells (Fig. 7). Differences in antigen presentation kinetics of peptides from the same viral protein have been shown  to account for the magnitude of expansion of antigen-specific CD8 T cells 56 , while epitope abundance has not 57 . Further studies are required to investigate these mechanisms in the pig IAV infection model, and the importance of the different epitopes for protective immunity.
Mechanisms by which TRM populations persist have been investigated by several authors in mice with sometimes contradictory results, so that it remains unclear whether some or all respiratory tract TRM populations are maintained by recruitment of circulating cells. While some experiments suggest that antigen encounter in the lung environment is important for TRM establishment and maintenance 42,53,58,59 , it has been reported that lung TRM may develop independently of pulmonary antigen encounter under specific inflammatory conditions 60 . Our results indicate that TEM but not TCM can be recruited to the lungs and since relative numbers of TEM in the circulation decrease with time after antigen exposure while TCM increase, (Fig. 4b) this may explain why recruitment of TRM declines with time 60 . These data are in accord with parabiosis experiments in mice since if the parabiosis is carried out early after immunization, (activated) antigen-specific cells from the immunized animal enter the lungs of the parabiotic partner but not if the parabiosis is performed later 59 .
Our data indicates that following influenza infections an airway population of TRM, recoverable by BAL, is established with a peak at 20-30 DPI of predominantly AAV-specific cells. We hypothesize that the late dominance of AAV may be because cells with this specificity continue to divide for longer than those specific for DFE or VAY, with continuing generation of cells with a TEM phenotype that can enter the lungs. Early on, BAL TRM cells express genes characteristic of T cell activation and 5-22% express Ki67, though Ki67 expression is much more prominent in tetramer binding cells in TBLN and PBMC, suggesting that the bulk of the BAL population arises by cell division prior to entry into the airways. BAL Ki67 expression declines over time and, since the tetramer binding population also declines, cell division clearly does not maintain the population, nor can it be replenished from the blood since the frequency of tetramer binding cells in blood declines similarly. The low level of Ki67 expression is in line with recent studies showing lower frequencies of Ki67 + TRM in the lung of healthy human donors and in SARS-CoV-2 infection compared to circulating T cells 16,61,62 .
Interestingly, expression of the transcription factors T-bet and Eomes has similar kinetics to Ki67 in BAL. In mice, T-box transcription factor downregulation and the consequent expression of CD103 on CD8 T cells are essential for TRM formation 63 . Although the lack of phenotypic changes over time may suggest that airway TRM are a stable effector population, this is not the case since there are major change in gene expression, with genes involved in cytokine production, differentiation and the TGFβ pathway highly upregulated at 21 compared to 7 DPI and expression declining at 63 DPI. TGFβ plays a critical role in tissue retention and differentiation into TRM [64][65][66] suggesting that tetramer binding cells at 21 DPI already present TRM features. In addition, we identified different sets of genes characteristic of TRM in humans, which are also upregulated at each of these time points. Viral antigen has been shown to persist in local lymph nodes for more than 2 months after IAV infection and in the lung for 110 days after Adenovirus vectored immunization in mice, influencing the recruitment of circulating antigen specific T cells to the lung and TRM gene expression 67,68 . In addition, others have reported marked differences in gene expression between CD8 T cells specific for different epitopes 69 , which should be investigated in the future to understand the prevalence of AAV + CD8 T cells.
These data provide a detailed analysis of the phenotype and distribution of tetramer binding CD8 T cells following IAV infection but not the protective efficacy of these cells in the pig model. However, when S-FLU immunized pigs are challenged with a heterosubtypic IAV strain, in the absence of cross-neutralizing antibodies, lung pathology is reduced but the viral load and viral shedding are not, in contrast to results in ferrets 27 and mice 70,71 . Surprisingly, after influenza infection the frequency of CD8 tetramer binding cells declines rapidly so that by day 63 there are very low frequencies in all tissues we have examined, except BAL. Our study was limited to CD8 T cells but future work will need to establish if other memory populations decline with similar kinetics and whether this decline correlates with loss of protection against lung pathology following heterosubtypic challenge.
Nevertheless, it appears that T cell memory to IAV may decline rapidly in this natural host species not maintained under specific pathogen free conditions, contrasting with results in transgenic mouse models 10 . Whether the pig or small animals better represent the way in which humans respond to influenza viruses remains to be determined. The fact that humans are susceptible to new influenza strains does suggest that heterosubtypic immunity declines with time and is not effective in preventing infection or transmission, although perhaps providing some degree of protection against severe disease 72 .The similarities in phenotype and transcriptional profile of porcine and human TRM highlight the value of this large animal model for understanding the importance of specificity of response in establishing protection against respiratory infections and for the development of next generation vaccines.

MATERIALS AND METHODS Animals and influenza challenge experiment
Animal experiments were conducted according to the UK Government Animal (Scientific Procedures) Act 1986 at the University of Bristol, the Pirbright Institute (TPI) and the Animal and Plant Health Agency (APHA) under project license P47CE0FF2, with approval from ethics committees at each institute. All institutions conform to the ARRIVE guidelines. The inbred Babraham pigs 33,73 were bred at APHA in line with the principles outlined in the FELASA recommendations 74 . Animals were pre-screened to ensure absence of exposure to influenza viruses by hemagglutination inhibition test. The pigs (female and male) were randomized and acclimatized for at least 7 days before any procedures were carried out.
Three influenza challenge experiments were performed at the University of Bristol (T1, T2 and T3), as previously reported 24 while a longer time course study (T4) took place at APHA (Fig. 2a). For T1, T2 and T3 thirty-eight Babraham inbred pigs (9.3 weeks old on average) were experimentally infected intranasally, using a mucosal atomisation device (MAD) (MAD300, Wolfe-Tory Medical), with 1 × 10 7 PFU of MDCK grown H1N1 A/swine/ England/1353/2009 (H1N1pdm09), 2 ml per nostril. In T1 and T2 one infected pig was culled each day on days 1 to 7, 9, 11 and 13 post challenge while two control uninfected pigs were sampled prior to infection and two more at day 8. During T3 three challenged animals were culled at days 6, 7, 13, 14, 20 and 21 post infection. Six control uninfected animals were included in T3: three culled at day −1 and the other three on the day of challenge. During T4 twelve pigs (10.2 weeks old on average) were challenged as described above and four pigs were randomly culled at each of day 21, 42 and 63 post infection. All animals were included in the analysis and no criteria for exclusion were established a priori. A first immunization experiment was performed at TPI where three Babraham pigs of 5.5 weeks of age were sedated and immunized with H7N1 S-FLU [eGFP/N1(PR8)] H7t(Netherlands/219/2003) (2.4. x 10 8 50% tissue culture infective dose (TCID50)) administered by aerosol using a SOLO vibrating mesh nebuliser (Aerogen ltd.) as previously described 25 . A group of two control pigs was left untreated. Vaccinated pigs were boosted after 3 weeks as described above and culled 3 weeks post boost (WPB). In total, 10 min prior to sacrifice animals were infused intravenously (i.v.) with anti-CD3 monoclonal antibody (mAb) (clone PPT3), produced in house, at a concentration of 1 mg/kg to label circulating T cells. A second immunization experiment was performed at APHA. Six Babraham pigs received H1N1 S-FLU [eGFP*/N1(A/Eng/195/2009)] H1(A/Eng/195/2009) (7.3 × 10 7 TCID50) by aerosol twice, three weeks apart, as described before 25 . Five unimmunized pigs were used as controls. One pig in the control group reached its humane end point because of a pre-existing heart condition, limiting the number of pigs in the control to five animals. At 3WPB, all animals were challenged intranasally with 2.8 × 10 6 PFU of H1N1pdm09 delivered by MAD and culled four days later. To minimise potential confounders pigs were challenged with H1N1pdm09 at the same time, in random order not based on prior treatment.
Half of the pigs in each group received anti-CD3 mAb i.v. 10 min prior to culling.
The time point selected (4 days post challenge) is based on our extensive experience and published work with other swine/zoonotic influenza viruses to provide the best snap shot of disease development.
As each group was culled at different time or immunized or not, it was not possible to run completely blinded experiments.

IFNγ ELISpot assay and identification of influenza nucleoprotein epitopes
Cryopreserved lymphocytes were stimulated using H1N1pmd09 (MOI = 1), medium and H1N1pmd09 nucleoprotein derived peptides (GL Biochem Ltd.) and frequencies of IFNγ spot forming cells were determined as described before 75 . To identify CD8 T cell epitopes in H1N1pdm09 nucleoprotein pools of ten 18 amino acid (aa) peptides overlapping by 12 aa, were used for initial ELISpot screening. Individual peptides from responding pool were then used to identify the highest responding peptides within each pool. Minimal epitopes were then defined using a second screen with 9 aa peptides derived from the 18 aa peptides giving the highest responses ( Supplementary Fig. 2). Pool 3 and 4 consistently showed the highest responses across tissues and therefore were broken down to identify minimal epitopes of 9 aa length: NP 181-189 AAVKGVGTI (AAV) and NP 217-225 VAYERMCNI (VAY), which were confirmed to be CD8 epitopes (data not shown) and loaded into porcine SLA-2 molecules to generate tetramers. Responses to pool 5 were attributed to the previously defined DFE epitope and this pool was therefore not further deconvoluted 34 .

Flow cytometric analysis and cell sorting
Cryopreserved single cell suspensions from PBMC, TBLN, lung, BAL, NT and trachea were thawed and rested for 2 hours in RPMI supplemented with Glutamax, 1% Penicillin-Streptomycin, 5% HEPES and 10% FBS (all from Thermo Fisher) at room temperature before staining. Two million cells per well were stained in 96 well plates with each NP tetramer (AAV, DFE (NP 290-298 DFEREGYSL ) and VAY) separately (as previously described 25,34 ). Following tetramer staining, fluorochrome-conjugated antibodies (see Supplementary Table 5  CD45RA -CCR7 -) using a BD FACSAria cell sorter. Cells were then centrifuged at 1000 × g for 5 min and re-suspended in cell culture medium (RPMI, 10% FBS, HEPES, Sodium pyruvate, Glutamax and Penicillin/Streptomycin) overnight at 37°C. On the following day, cells were stimulated using PMA Ionomycin (BioLegend) for 2, 4 and 6 hours and cytokines detected by intracellular cytokine staining (see Supplementary Table 5 for antibodies list) using Fixation/Permeabilization Solution Kit (BD Biosciences), unstimulated cells were used as a control. Cells isolated from BAL, TBLN and PBMC were stimulated with H1N1pdm09 (MOI = 1) overnight and cytokine production quantified by intracellular cytokine staining as previously described 24 . Data were analyzed by Boolean gating. Media only controls for each sample were used as baseline and subtracted from the results.
For preparation of RNA, BAL cells isolated from 7, 21 and 63 DPI samples (n = 3) were depleted of alveolar macrophages using a MACS cell separation LS column (Miltenyi Biotec) after staining for CD14 and CD172a (see Supplementary Table 5). Lymphocytes were then stained for DFE tetramer, CD8 and live/dead marker using a BD FACSAria (BD Bioscience).

RNA-seq
RNA-sequencing was performed on DFE + T cells from BAL because the tissue samples were sorted before the AAV tetramer was produced. BAL DFE + T cells were sorted (average of 4480 cells/sample) in PBS, samples were then centrifuged at 3000 g for 5 min. RNA was extracted using PicoPure RNA Isolation Kit (Thermo Fisher) according to the manufacturer instructions followed by DNase treatment (TURBO DNA-free™ Kit, Thermo Fisher). Isolated RNA was used as input for SMARTer Stranded Total RNA-Seq Kit v3 -Pico Input Mammalian (Takara Bio) and PCR amplification performed. cDNA was pooled and sequenced on NovaSeq using an S1 100 PE flow cell. Raw fastq files were used for an initial quality control using FastQC (https://www. bioinformatics.babraham.ac.uk/projects/fastqc/). Next, CogentAP (Cogent NGS Analysis Pipeline v1.0, Takara Bio) was used to trim and add the sample barcodes to the fastq header for each sample. The reads were then trimmed of Illumina and library prep adapters using cutadapt 77 , and subsequentially aligned to Sus scrofa Sscrofa11.1 assembly (GCA_000003025.6) using STAR 78 with default parameters. UMI-tools was used to discard duplicated reads 79 and featureCounts extracted the number of reads aligned to each gene feature 80 . Differential gene expression (DGE) analysis was carried out using TCC-GUI, which iterates DESeq2 for data normalisation, using the feature-Counts output for pairwise comparisons 81,82 . Pathway analysis was performed using Webgestalt online tool 83 , selecting GSEA as comparison method and sscrofa as reference genome. KEGG pathway database was used as reference and ensemble ID were uploaded together with ranked score (-log 10 (p value)*log 2 (fold change difference)) based on the results of DGE comparison. Other settings include: minimum number of IDs in the category (5), maximum number of IDs in the category (2000), significance level (FDR < 0.05) and number of permutation (1000).

Statistical analysis
To compare the changes over time in the proportion of CD8 + T cells specific for each tetramer in a tissue, the following curve was fitted to the data, namely, where y is the proportion at t days post infection, y 0 is the baseline proportion, y max is the peak proportion and t max is the time of peak proportion. This curve was chosen as it gives an appropriate shape for the data without including a large number of parameters. Note that the decay rate from peak titre is given by 6/t max . Similarly, changes over time in the proportion of tetramer-specific T cells expressing different cell surface markers in each tissue were compared by fitting exponential curves to the data, yðtÞ ¼ y 0 expðÀdtÞ; where y is the proportion at t days post infection, y0 is the initial proportion and d is the decay rate (/day). Each marker was analyzed independently.
In both analysis, variation in parameters (i.e. y 0 , y max and t max or y 0 and d) amongst tissues and tetramers was assessed by fitting different models to the data by nonlinear least squares and comparing the residual deviance for the models using F-tests 84 (Supplementary Tables 1 and 4). These analyses were implemented in Matlab (version R2020b; The Mathworks, Inc.).
Because data on CD8 + T cells in the trachea were only available for a reduced number of time points, the proportion of cells specific for each tetramer in this tissue were compared at each time point using a Kruskal-Wallis test followed by pairwise Wilcoxon rank-sum tests. This analysis was implemented in R (version 4.0.5) 85 .
To assess changes over time in the relative proportion of CD8 + T cells specific for each tetramer, the following model was used: where y(t) is the proportion specific for the tetramer t days post infection, K1 and K2 are the minimum and maximum frequencies, b the rate of change in frequency and d is the time of the maximum rate of change. This formulation ensures the total proportion of cells is 100%. Variation in parameters (i.e. K 1 , K 2 , b and d) amongst tissues was assessed by fitting different models to the data by nonlinear least squares and comparing the residual deviance for the models using F-tests 84 (Supplementary Tables 2). This analysis was implemented in Matlab (version R2020b; The Mathworks, Inc.). Only samples for which the frequency of all three tetramers was available were included in this analysis.
Trends in the phenotype of tetramer specific CD8 + T cells in each tissue were assessed using linear models. Each model included the proportion of cells in the population (CD45RA + / − and CCR7 + / − ) specific for a tetramer as the response variable and tissue, tetramer and days post infection as explanatory variables, together with two-and three-way interactions between the explanatory variables. Model simplification proceeded by stepwise deletion of non-significant (P > 0.05) terms (as judged by F-tests). This analysis was implemented in R (version 4.0.5) (R Core Team 2021).
Frequencies of cytokine secreting cells in DFE + and AAV + cells over time were compared using two-way ANOVA in GraphPad Prism version 9.1.0.
The number of pigs in each experiment are indicated in the figure legends. Sample size for the immunization experiment was calculated based on previous experiments, with a SD of approximately 0.5 and 80% power and 95% confidence 27,75 .